En 2018, le PSAR analyse urbaine, ancêtre de la section analyse urbaine (DG/DDAR/DAR/DSAU), a développé un package R, nommé btb.
Sa principale fonction, kernelSmoothing, permet de réaliser très facilement un carroyage et un lissage sur des données géolocalisées. Grâce à différentes extensions, ce package peut être utilisé de plusieurs façons :
BTBpyÀ partir de données ponctuelles, nous allons apprendre :
Les traitements sont réalisés avec le langage R. Le code de cette présentation est disponible ici.
Avant toute diffusion auprès des partenaires, il faut bien veiller à respecter :
Pour simplifier les programmes présentés dans ce TP, les représentations graphiques ne respectent pas les règles élémentaires de la sémiologie cartographique très bien résumées par cette infographie :
Auteur : Timothée Giraud, auteur de la librairie mapsf
Il faudra bien veiller à appliquer ces règles de sémiologie sur les cartes que vous réaliserez ultérieurement.
Pour information, voici les systèmes de projection que vous pouvez régulièrement rencontrer pour la métropole :
| Nom | Description | Code EPSG |
|---|---|---|
| Lambert93 | Système de projection officiel pour la métropole | 2154 |
| LAEA | Système de projection européen | 3035 |
| WGS84 | GPS (utile pour utiliser Leaflet) | 4326 |
Quatre librairies sont nécessaires pour ce TP.
sf pour manipuler des fichiers spatiaux (importer des .shp, transformer des projections, et réaliser des géotraitements)mapsf pour réaliser des cartes dans RStudioleaflet pour réaliser des cartes interactives (fond de carte OpenStreetMap)btb pour le carroyage et lissagedplyr pour le traitement des données, en particulier l’agrégation géographique.Remarque 2 : Le choix de dplyr plutôt que data.table se justifie ici du fait de sa forte compatibilité avec les objets géomatiques.
# Charger les librairies nécessaires
## Liste des librairies utilisées
packages <- c("dplyr","sf","btb","mapsf","leaflet")
## Vérifier si la librairie est installée, si non l'installer, puis la charger
package.check <- lapply(
packages,
FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)
}
}
)
Pour ce TP, nous utilisons la base « Demandes de Valeurs Foncières », ou DVF, qui recense l’ensemble des ventes de biens fonciers réalisées au cours des cinq dernières années, en métropole et dans les départements et territoires d’outre-mer — sauf à Mayotte et en Alsace-Moselle. Les biens concernés peuvent être bâtis (appartement et maison) ou non bâtis (parcelles et exploitations). Les données sont produites par la direction générale des finances publiques. Elles proviennent des actes enregistrés chez les notaires et des informations contenues dans le cadastre.
A partir de cette source, nous avons constitué une base de données qui s’intéresse uniquement au périmètre de la petite couronne parisienne (départements 75, 92, 93 et 94) et au millésime 2021.
La base contient les 8 variables suivantes :
id_mutation : identifiant unique de la ventedate_mutation : date de la ventetype_local : appartement ou maisonnombre_pieces_principales : Nombre de pièces dans le logementvaleur_fonciere : prix de ventesurface_reelle_bati : surface du logementx : longitude (en projection Lambert 93)y : latitude (en projection Lambert 93)Le code ayant permis de constituer la base de données est disponible ici : XXXXXX. TODO. + préciser éventuellement à grands traits les filtres réalisés dans la remarque ci-dessous.
Remarque importante : Cette base a été filtrée de manière à être la plus pédagogique possible pour cette formation. Elle n’est donc ni représentative de la réalité ni fidèle au champ proposé par le producteur. Mieux vaut donc de pas appuyer vos investissements immobiliers sur les résultats de ce TP.
Le code ci-dessous permet d’importer la première base de données ventesImmo_couronneParis.RDS utilisée dans ce tutoriel. Elle est stockée sous Minio, dans le “bucket public” : s3/kantunez/diffusion/projet_formation/r_lissage_spatial/.
# Charger la source de données (variable nommée dfBase) avec S3
url_bucket <- "https://minio.lab.sspcloud.fr/"
bucket <- "kantunez"
object <- "diffusion/projet_formation/r_lissage_spatial/ventesImmo_couronneParis.RDS"
dfBase <-
aws.s3::s3read_using(
FUN = base::readRDS,
object = object,
bucket = bucket
,
opts = list("region" = "")
)
Il est aussi possible de charger cette base en dehors du SSPCloud en la téléchargeant grâce à son URL puis en la chargeant dans R.
#Se positionner sur le répertoire de travail.
#Non nécessaire si projet RStudio créé au préalable
setwd("Mon/repertoire/qui/va/bien")
# Charger la source de données (variable nommée dfBase) avec son URL
download.file(paste0(url_bucket,bucket,"/",object), destfile = "ventesImmo_couronneParis.RDS")
dfBase <- readRDS("ventesImmo_couronneParis.RDS")
On peut ensuite manipuler notre base de données chargée.
# Visualiser les premières lignes de la base
head(dfBase)
## id_mutation date_mutation type_local nombre_pieces_principales valeur_fonciere
## 1 2021-447023 2021-01-08 Appartement 3 480000
## 2 2021-447024 2021-01-05 Appartement 2 345000
## 3 2021-447025 2021-01-08 Appartement 3 384000
## 4 2021-447027 2021-01-08 Appartement 3 261900
## 5 2021-447029 2021-01-05 Appartement 2 407200
## 6 2021-447030 2021-01-15 Appartement 3 1040000
## surface_reelle_bati x y
## 1 64 647357.3 6868635
## 2 43 644483.5 6867695
## 3 41 648001.8 6866153
## 4 44 647414.8 6868937
## 5 24 646929.9 6864730
## 6 90 645132.7 6864646
dim(dfBase)
## [1] 34489 8
Nous allons charger le contour géographique des départements de la petite couronne parisienne. En cartographie, ces fichiers sont appelés des « couches vectorielles » (à ne pas confondre avec les couches « raster » qui correspondent, pas exemple, aux fonds de carte satellite OpenStreetMap ou Google Maps). Une fois les départements chargés, nous allons sélectionner le contour de la commune de Paris.
Différents types de fichiers permettent de stocker des couches vectorielles. Les deux principaux sont le Geopackage (.gpkg) et le Shapefile (.shp). Nous recommandons l’utilisation du .gpkg pour les raisons suivantes :
Pour aller plus loin sur la comparaison des formats, voir ici.
Récupérons le contour du territoire étudié avec la fonction read_sf du package sf.
# Charger le fond de carte du territoire étudié
object <- "diffusion/projet_formation/r_lissage_spatial/depCouronne.gpkg"
depCouronne_sf <- st_read(paste0(url_bucket,bucket,"/",object))
## Reading layer `depCouronne_sf' from data source
## `https://minio.lab.sspcloud.fr/kantunez/diffusion/projet_formation/r_lissage_spatial/depCouronne.gpkg'
## using driver `GPKG'
## Simple feature collection with 4 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
## Projected CRS: RGF93 / Lambert-93
# Visualisons cette couche vectorielle
head(depCouronne_sf)
## Simple feature collection with 4 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
## Projected CRS: RGF93 / Lambert-93
## code libelle reg surf geom
## 1 75 Paris 11 105 MULTIPOLYGON (((660897 6860...
## 2 92 Hauts-de-Seine 11 176 MULTIPOLYGON (((648796 6847...
## 3 93 Seine-Saint-Denis 11 237 MULTIPOLYGON (((659428 6861...
## 4 94 Val-de-Marne 11 245 MULTIPOLYGON (((656908 6846...
plot(depCouronne_sf$geom)
depCouronne_sf <- depCouronne_sf %>% rename(geometry=geom)
# Sélection de Paris
paris_sf <- depCouronne_sf[depCouronne_sf$code=="75",]
# Visualisons cette couche vectorielle
head(paris_sf)
## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
## Projected CRS: RGF93 / Lambert-93
## code libelle reg surf geometry
## 1 75 Paris 11 105 MULTIPOLYGON (((660897 6860...
plot(paris_sf$geometry)
Pour être certain que les données et le territoire soient dans le même système de projection, il est possible de transformer ce dernier à l’aide de la fonction st_transform.
# Transformer la projection
depCouronne_sf <- sf::st_transform(depCouronne_sf, crs = 2154)
paris_sf <- sf::st_transform(paris_sf, crs = 2154)
À partir de ces territoires, nous pouvons déterminer le rectangle englobant pour réduire la quantité de données manipulées.
Pour éviter les effets de bord, il faut également sélectionner des données au-delà de notre zone d’intérêt (ici Paris intramuros). Autrement dit, si on ne sélectionne que les ventes immobilières situées dans Paris intramuros, et qu’on lisse ce nuage de points, les zones situées à proximité du périphérique seront artificiellement peu denses en ventes immobilières. En effet, elles sont situées à proximité de zones artificiellement vides en ventes immobilières dans notre nuage de points, à savoir dans les communes limitrophes. Pourtant, dans la réalité, il y a bien des ventes réalisées à Malakoff, Vanves, Montrouge, etc. Et ces ventes influencent la densité lissée aux niveau de la Porte de Vanves.
# Mise en forme de la couche "territoire" : sélection des variables et renommage
paris_sf$nom <- "territoire"
paris_sf <- paris_sf[,c("nom","geometry")]
# Création d'une bbox vectorielle autour du territoire
bbox <- sf::st_bbox(paris_sf)
bbox
## xmin ymin xmax ymax
## 643076 6857499 660897 6867034
bbox_sf = st_sf(nom="bbox", geometry = st_as_sfc(bbox), crs = 2154)
bbox_sf
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
## Projected CRS: RGF93 / Lambert-93
## nom geometry
## 1 bbox POLYGON ((643076 6857499, 6...
# Création d'un buffer de la bbox, avec une marge de 2000 mètres
marge <- 2000
bufferBbox <- bbox
bufferBbox[["xmin"]] <- bufferBbox[["xmin"]]-marge
bufferBbox[["xmax"]] <- bufferBbox[["xmax"]]+marge
bufferBbox[["ymin"]] <- bufferBbox[["ymin"]]-marge
bufferBbox[["ymax"]] <- bufferBbox[["ymax"]]+marge
bufferBbox_sf = st_sf(nom="buffer_bbox", geometry = st_as_sfc(bufferBbox), crs = 2154)
# Assemblage des 3 couches dans une même table sf pour la cartographie
couches_plot <- rbind(bbox_sf,
bufferBbox_sf,
paris_sf)
# Passage d'une couche de polygones à une couche de lignes (pour la cartograhie)
couches_plot <- st_cast(couches_plot,"MULTILINESTRING")
# Cartographie pédagogique avec le package mapsf
mf_init(x = couches_plot, theme = "agolalight")
mf_map(x = st_cast(couches_plot[couches_plot$nom=="territoire",],"POLYGON"),
type = "base",
col="wheat",
border=FALSE,
add = TRUE)
mf_map(x = couches_plot,
type = "typo",
var = "nom",
pal = c("aquamarine4", "yellow3","wheat"),
lwd=8,
leg_title = "Zones d'analyse",
add = TRUE)
mf_layout(title = "Gestion du territoire d'analyse",credits = "Insee-DSAU, IGN, mapsf")
Ce schéma présente :
Pour sélectionner les données individuelles d’intérêt pour le lissage, on peut tout d’abord filtrer les observations dont les coordonnées sont comprises à l’intérieur du rectangle jaune. Ce filtre est très efficace computationnellement car il ne dépend que des valeurs numériques prises par les variables de longitude et de latitude et n’utilise pas les propriétés vectorielles des données qui sont des attributs plus chronophages à utiliser.
# Repérer les coordonnées extrêmes, avec une marge (ici 2km)
bufferBbox
## xmin ymin xmax ymax
## 641076 6855499 662897 6869034
xMin <- bufferBbox["xmin"]
xMax <- bufferBbox["xmax"]
yMin <- bufferBbox["ymin"]
yMax <- bufferBbox["ymax"]
# Ne garder que les données dans le rectangle englobant le territoire étudié, sans traitement vectoriel !
dfBase_filtre <- dfBase[dfBase$x >= xMin & dfBase$x <= xMax & dfBase$y >= yMin & dfBase$y <= yMax, ]
dim(dfBase_filtre)
## [1] 24419 8
Il est possible de procéder autrement. On peut utiliser nos données individuelles comme un ensemble de points (approche vectorielle) et procéder à des intersections géographiques. Pour cela :
sf vectorielCette méthode, bien qu’élégante, peut-être lourdre d’un point de vu computationnel. Avec des données volumineuses, il convient au minimum de faire un premier filtrage avec la précédente méthode.
# Transformer le tableau de données filtré en objet géographique
sfBase <- sf::st_as_sf(dfBase, coords = c("x", "y"), crs = 2154)
# Création d'un buffer autour du territoire
buffer_sf <- st_buffer(paris_sf, dist = 2000)
# Repérer les indices des observations contenues dans notre buffer d'intérêt
indiceObsContenues <- unlist(sf::st_contains(buffer_sf, sfBase))
# Réduire la base aux seules observations dans le territoire
dfBase_filtre <- dfBase[indiceObsContenues, ]
# Mise en forme de la couche buffer
buffer_sf$nom <- "buffer"
# Assemblage des 2 couches dans une même table sf pour la cartographie
couches_plot <- rbind(buffer_sf,paris_sf)
# Passage d'une couche de polygones à une couche de lignes (pour la cartographie)
couches_plot <- st_cast(couches_plot, "MULTILINESTRING")
# Cartographie pédagogique avec le package mapsf
mf_init(x = couches_plot, theme = "agolalight")
mf_map(x = couches_plot,
type = "typo",
var = "nom",
pal=c("yellow3","wheat"),
lwd=8,
leg_title = "Zones d'analyse",
add = TRUE)
mf_map(x = st_cast(couches_plot[couches_plot$nom=="territoire",],"POLYGON"),
type = "base",
col="wheat",
border=FALSE,
add = TRUE)
# Échantillon de 500 observations dans le buffer
sfBase_sample <- sfBase[indiceObsContenues, ]
sfBase_sample <- sfBase_sample[sample(1:nrow(sfBase_sample),2000) ,]
mf_map(sfBase_sample,
add = TRUE,
leg_title = "observations")
mf_layout(title = "Gestion du territoire d'analyse",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Remarque: Pour la zone tampon, il convient de prendre une marge légèrement plus grande que le rayon de lissage envisagé, ceci afin d’éviter les effets de bord tout en limitant les temps de calcul. Une fois le lissage réalisé, on peut réduire la carte à la zone étudiée.
Avant de lisser les données ponctuelles, on peut souhaiter représenter ces données sous forme carroyée afin de se les approprier. Il faut pour cela :
Le code ci-dessous procède de la sorte :
dfToGrid du package btb.La fonction dfToGrid attend comme paramètres obligatoires :
df : un tableau avec les colonnes x et y représentant les coordonnées des centroïdes de la grille. (Remarque : la grille n’est pas nécessairement régulière.)sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé.# Arrondir les coordonnées des observations aux centroides les plus proches
iCellSize = 200 # Carreaux de 200m
points_carroyage <- dfBase_filtre
points_carroyage$x_centroide = as.integer(floor(points_carroyage$x / iCellSize) * iCellSize + (iCellSize / 2))
points_carroyage$y_centroide = as.integer(floor(points_carroyage$y / iCellSize) * iCellSize + (iCellSize / 2))
head(points_carroyage)
## id_mutation date_mutation type_local nombre_pieces_principales valeur_fonciere
## 3 2021-447025 2021-01-08 Appartement 3 384000
## 5 2021-447029 2021-01-05 Appartement 2 407200
## 6 2021-447030 2021-01-15 Appartement 3 1040000
## 9 2021-447034 2021-01-18 Appartement 4 625000
## 10 2021-447035 2021-01-18 Appartement 4 721250
## 11 2021-447037 2021-01-07 Appartement 4 1035000
## surface_reelle_bati x y x_centroide y_centroide
## 3 41 648001.8 6866153 648100 6866100
## 5 24 646929.9 6864730 646900 6864700
## 6 90 645132.7 6864646 645100 6864700
## 9 68 647494.6 6866019 647500 6866100
## 10 94 643891.9 6864587 643900 6864500
## 11 94 648239.0 6866235 648300 6866300
# Compter le nombre de ventes par carreau (si on cherche à représenter la densité)
points_carroyage <- points_carroyage %>%
group_by(x=x_centroide,y=y_centroide) %>%
count(name = "nbVentes")
# Générer la grille
points_carroyage <- btb::dfToGrid(df = points_carroyage, sEPSG = "2154", iCellSize = iCellSize)
# Restriction du champ : on ne retient que les carreaux intersectant Paris
points_carroyage <- points_carroyage[unlist(st_intersects(paris_sf,points_carroyage)),]
# Cartographie
mf_init(x=points_carroyage,theme = "agolalight")
mf_map(x = points_carroyage,
type = "choro",
var="nbVentes",
breaks = "quantile",
nbreaks = 5,
lwd=1,
leg_val_rnd = 1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Carroyage du nombre de ventes",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Remarque 1 : Le carroyage peut aussi être utilisé pour simplifier les données avant lissage si le calcul de lissage est trop long en raison d’un très grand nombre d’observations.
Remarque 2 : Les étapes avant la génération de la grille ne sont pas intégrées dans des fonctions du package btb, il faut donc s’inspirer du code ci-dessus pour reproduire le même type de carroyage. Il est envisagé qu’une future version du package btb simplifie le code permettant de réaliser un carroyage.
Un lissage simple à réaliser correspond au calcul d’une densité, à savoir une quantité par unité de surface (un carreau). Dans cette section nous allons calculer la densité ventes de logements sur la ville de Paris au cours de l’année 2021.
Pour ce faire, nous allons au préalable créer une variable nbObsLisse qui vaudra 1 pour chaque observation (donc logement vendu en 2021) et qui sera lissée avec la fonction kernelSmoothing du package btb.
Toujours pour éviter les effets de bord, le lissage sera effectué sur une zone plus large que la zone d’intérêt (utilisation d’un buffer).
Pour ce premier lissage, nous choisissons des carreaux de 200 mètres et un rayon de lissage de 400 mètres.
La fonction de lissage btb::kernelSmoothing permet de générer automatiquement le carroyage préalable au lissage.
# Créer une nouvelle variable pour compter le nombre de ménages lissé
dfBase_filtre$nbObsLisse <- 1
# Visualiser les premières lignes de la base
head(dfBase)
## id_mutation date_mutation type_local nombre_pieces_principales valeur_fonciere
## 1 2021-447023 2021-01-08 Appartement 3 480000
## 2 2021-447024 2021-01-05 Appartement 2 345000
## 3 2021-447025 2021-01-08 Appartement 3 384000
## 4 2021-447027 2021-01-08 Appartement 3 261900
## 5 2021-447029 2021-01-05 Appartement 2 407200
## 6 2021-447030 2021-01-15 Appartement 3 1040000
## surface_reelle_bati x y
## 1 64 647357.3 6868635
## 2 43 644483.5 6867695
## 3 41 648001.8 6866153
## 4 44 647414.8 6868937
## 5 24 646929.9 6864730
## 6 90 645132.7 6864646
Pour réaliser le lissage, il suffit de fournir à la fonction kernelSmoothing quatre paramètres :
dfObservations : le tableau des données à lisser. Il doit nécessairement contenir une colonne x, une colonne y, et 1 à n colonnes numériques (variables à lisser) ;sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé ;iCellSize : un entier indiquant la taille des carreaux ;iBandwidth : un entier indiquant le rayon de lissage.Attention, la fonction retourne une erreur en cas de :
# Lisser
dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 400)
# Visualiser les premières lignes des données lissées
head(dfBaseLisse)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 650200 ymin: 6855400 xmax: 653400 ymax: 6855800
## Projected CRS: RGF93 / Lambert-93
## x y nbObsLisse geometry
## 1 651500 6855500 1.0477469 POLYGON ((651400 6855600, 6...
## 2 652100 6855500 0.8521979 POLYGON ((652000 6855600, 6...
## 3 652700 6855500 2.7218839 POLYGON ((652600 6855600, 6...
## 4 653300 6855500 1.4556305 POLYGON ((653200 6855600, 6...
## 5 650300 6855700 2.1375004 POLYGON ((650200 6855800, 6...
## 6 650700 6855700 2.0552302 POLYGON ((650600 6855800, 6...
# Nombre de lignes lissées
nrow(dfBaseLisse)
## [1] 3259
La fonction kernelSmoothing génère une grille de carreaux, et chaque carreau est affecté de la densité lissée en son centroïde. À noter que la grille carroyée épouse globalement le périmètre des points en entrée de la fonction (elle va même un peu au-delà par défaut, cf. séquence théorique sur le lissage de la formation “Comment utiliser les outils de l’analyse urbaine ?”).
mf_init(x=paris_sf,theme = "agolalight")
mf_map(x = st_cast(dfBaseLisse[,c("geometry")],"LINESTRING"),
lwd=1,add=T)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=8,
col="wheat",add = TRUE)
mf_layout(title = "Gestion du territoire d'analyse",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Une propriété particulièrement importante de la fonction de lissage classique du package btb est qu’elle est conservative. Cela signifie que nous avons la même somme des variables additives sur le champ géographique concerné avant et après lissage.
sum(dfBase_filtre$nbObsLisse) # Nombre de ventes dans la petite couronne
## [1] 22221
sum(dfBaseLisse$nbObsLisse) # Après lissage, nombre lissé de ventes dans les carreaux
## [1] 22221
Remarque : Veillez toujours à respecter le secret statistique au moment de la publication de vos cartes ! Par exemple, si vous utilisez des données issues des bases fiscales, vous ne pouvez représenter des informations sur des carreaux comportant moins de 11 observations, même si ce nombre est lissé… Pour ce faire, vous pouvez réaliser le filtre ci-dessous :
# Exclure les carreaux ne respectant pas le secret statistique
dfBaseLisse <- dfBaseLisse[dfBaseLisse$nbObsLisse >= 11, ]
Voyons maintenant le nombre lissé de ventes de logements à Paris en 2021 suite à ce premier lissage.
Quelques informations préalables :
dfBase_filtre contient toutes les ventes comprises dans un buffer de 2 km autour de Paris.sf::st_intersects).choroplèthe ;# Filtrage des carreaux lissés intersectant la ville de Paris
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
# Carte lissée
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 400m",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Remarque : la variable lissée nbObsLiss correspond au nombre de ventes (observations) par carreau (ici de 200 mètres). Ainsi, il est possible d’obtenir une densité au km² en multipliant la variable nbObsLiss par 25 dans le cas présent.
Il n’est pas possible de connaitre a priori la taille optimale du rayon de lissage. Il faut donc en essayer plusieurs pour trouver le meilleur compromis entre précision et généralisation. Nous allons donc faire le même lissage que précédemment avec un rayon de 600 m au lieu de 400 m.
dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 600)
# Visualiser les premières lignes des données lissées
head(dfBaseLisse)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
## x y nbObsLisse geometry
## 1 651300 6855300 0.1065710 POLYGON ((651200 6855400, 6...
## 2 651500 6855300 0.1384200 POLYGON ((651400 6855400, 6...
## 3 651700 6855300 0.2001949 POLYGON ((651600 6855400, 6...
## 4 651900 6855300 0.2094593 POLYGON ((651800 6855400, 6...
## 5 652100 6855300 0.1441657 POLYGON ((652000 6855400, 6...
## 6 652300 6855300 0.1852248 POLYGON ((652200 6855400, 6...
# Filtrage des carreaux lissés dans Paris
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
# Carte lissée
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 600m",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
À noter qu’on peut réitérer l’opération avec un rayon de lissage de 1000 mètres.
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 651000 ymin: 6855000 xmax: 652200 ymax: 6855200
## Projected CRS: RGF93 / Lambert-93
## x y nbObsLisse geometry
## 1 651100 6855100 0.09890581 POLYGON ((651000 6855200, 6...
## 2 651300 6855100 0.09483879 POLYGON ((651200 6855200, 6...
## 3 651500 6855100 0.11358059 POLYGON ((651400 6855200, 6...
## 4 651700 6855100 0.14485863 POLYGON ((651600 6855200, 6...
## 5 651900 6855100 0.16973038 POLYGON ((651800 6855200, 6...
## 6 652100 6855100 0.18790912 POLYGON ((652000 6855200, 6...
À mesure que nous augmentons le rayon de lissage, la carte révèle des aspects structurels des données. Mais ceci se fait au détriment des spécificités locales visibles seulement avec un petit rayon de lissage. Le choix revient dont au statisticien-géographe, en fonction de sa connaissance des données et des objectifs recherchés.
Enfin, la grille carroyée est ici visible à des fins pédagogiques, il est bien-entendu possible et conseillé de la supprimer sur les cartes publiées comme dans cet exemple :
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="nbObsLisse",
breaks = "quantile",
nbreaks = 5,
border = NA, # C'est ici que ça se passe
# lwd=1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m, sans grille",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Ci-dessous, nous illustrons les effets de bord produits si on lisse uniquement les ventes réalisées dans Paris intramuros. Dans les deux cartes ci-dessous, le rayon de lissage (2 000 mètres) et les bornes de discrétisation sont communs pour assurer la comparabilité des deux cartes.
Ainsi, on peut remarquer que les effets de bord se matérialisent à la périphérie de Paris, où les densités lissées sont artificiellement plus faibles que quand on lisse en prenant une marge.
# Lissage avec gestion des effets de bord
dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 2000)
# Lissage sans gestion des effets de bord (Base des ventes uniquement dans intramuros)
sfBase_intramuros <- sfBase[unlist(st_contains(paris_sf,sfBase)),]
dfBase_intramuros <- sfBase_intramuros %>%
mutate(x=st_coordinates(geometry)[,1],
y=st_coordinates(geometry)[,2],
nbObsLisse=1) %>%
st_drop_geometry() %>%
select(nbObsLisse,x,y)
dfBaseLisse_intramuros <- btb::kernelSmoothing(dfObservations = dfBase_intramuros,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 2000)
Évoquer le paramétrage de la grille carroyée afin de prendre en compte les limites artificielles (limites administratives) et les limites physiques (mer, montagne, etc.)
Dans cette section nous allons nous intéresser au prix moyen des logements vendus à Paris en 2021.
Remarque : Une moyenne n’est pas le lissage du rapport (liss(variable / nbObsLisse)) mais le rapport des lissages : liss(variable) / liss(nbObs).
Ainsi, il convient de :
Dans l’exemple ci-dessous, on prend un rayon de lissage de 600 m :
dataLissage <- dfBase_filtre[,c("nbObsLisse","valeur_fonciere","x","y")] # Ajout de la variable "valeur_fonciere"
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 600)
# Visualiser les premières lignes des données lissées
head(dfBaseLisse)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
## x y nbObsLisse valeur_fonciere geometry
## 1 651300 6855300 0.1065710 38517.44 POLYGON ((651200 6855400, 6...
## 2 651500 6855300 0.1384200 47025.68 POLYGON ((651400 6855400, 6...
## 3 651700 6855300 0.2001949 66234.22 POLYGON ((651600 6855400, 6...
## 4 651900 6855300 0.2094593 72332.24 POLYGON ((651800 6855400, 6...
## 5 652100 6855300 0.1441657 52248.55 POLYGON ((652000 6855400, 6...
## 6 652300 6855300 0.1852248 68258.16 POLYGON ((652200 6855400, 6...
# Calculer le taux
dfBaseLisse$prixMoyen <- dfBaseLisse$valeur_fonciere / dfBaseLisse$nbObsLisse
# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="prixMoyen",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Lissage des prix de vente avec un rayon de 600m",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
On remarque des prix moyens particulièrement élevés dans le centre et l’Ouest de la capital, ce qui semble conforme à l’intuition.
Attention cependant aux éventuelles valeurs atypiques (ventes extrêmement chères, erreur dans les données) : le lissage en moyenne est potentiellement sensible à ces valeurs atypiques (contrairement au lissage quantile, voir infra).
Dans cette section, nous allons nous intéresser à la proportion de ventes portant sur des grands logements (plus de 4 pièces principales).
# Créer une nouvelle variable pour détecter les logements de plus de 4 pièces principales
dfBase_filtre$quatrePieces <- ifelse(dfBase_filtre$nombre_pieces_principales >= 4, 1, 0)
dataLissage <- dfBase_filtre[,c("nbObsLisse","quatrePieces","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 600)
# Visualiser les premières lignes des données lissées
head(dfBaseLisse)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
## x y nbObsLisse quatrePieces geometry
## 1 651300 6855300 0.1065710 0.07578049 POLYGON ((651200 6855400, 6...
## 2 651500 6855300 0.1384200 0.07838126 POLYGON ((651400 6855400, 6...
## 3 651700 6855300 0.2001949 0.05794216 POLYGON ((651600 6855400, 6...
## 4 651900 6855300 0.2094593 0.02832061 POLYGON ((651800 6855400, 6...
## 5 652100 6855300 0.1441657 0.02332406 POLYGON ((652000 6855400, 6...
## 6 652300 6855300 0.1852248 0.13111821 POLYGON ((652200 6855400, 6...
# Calculer le taux
dfBaseLisse$txQuatrePieces <- dfBaseLisse$quatrePieces / dfBaseLisse$nbObsLisse
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
# Cartographie
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="txQuatrePieces",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Logements avec 4 pièces ou plus (rayon de 600m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
La carte des prix moyens et la carte des grands logements se ressemblent, fort logiquement. On peut alors avoir envie de lisser le prix au mètre-carré.
Pour ce faire, n’oubliez pas qu’il convient de lisser séparément les prix de vente (numérateur) et le nombre de mètres-carrés (dénominateur) des logements vendus. Autrement, si on lisse le prix au m² de chaque logement vendu, on sur-pondère artificiellement les prix au m² des petits logements (car l’unité statistique devient alors le logement, et non le m² vendu).
dataLissage <- dfBase_filtre[,c("surface_reelle_bati","valeur_fonciere","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
sEPSG = "2154",
iCellSize = 200,
iBandwidth = 600)
# Visualiser les premières lignes des données lissées
head(dfBaseLisse)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
## x y surface_reelle_bati valeur_fonciere geometry
## 1 651300 6855300 8.067281 38517.44 POLYGON ((651200 6855400, 6...
## 2 651500 6855300 10.163131 47025.68 POLYGON ((651400 6855400, 6...
## 3 651700 6855300 13.385886 66234.22 POLYGON ((651600 6855400, 6...
## 4 651900 6855300 12.786188 72332.24 POLYGON ((651800 6855400, 6...
## 5 652100 6855300 7.943338 52248.55 POLYGON ((652000 6855400, 6...
## 6 652300 6855300 12.310115 68258.16 POLYGON ((652200 6855400, 6...
# Calculer le prix au m²
dfBaseLisse$prixM2 <- dfBaseLisse$valeur_fonciere / dfBaseLisse$surface_reelle_bati
# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="prixM2",
breaks = "quantile",
nbreaks = 5,
lwd=1,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Prix au m² (rayon de 600m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Dans cette section, nous allons calculer le 1er décile, la médiane et le 9ème décile du prix de vente des logements à Paris. Ce calcul se fait toujours avec la fonction kernelSmoothing mais il faut ajouter un paramètre, vQuantiles, qui contient la liste des quantiles souhaités.
Remarque : Le calcul de quantiles géographiquement pondérés présente certaines différences par rapport au lissage classique vu jusqu’à présent :
nbObs indiquant le nombre d’observations réel (non lissé) ayant contribué au calcul du carreau.On représente ici le lissage du premier décile et de la médiane.
Pour rappel, le lissage quantile (médiane, déciles . . . ) permet d’avoir des indicateurs moins sensibles aux valeurs extrêmes et d’enrichir l’analyse avec des indicateurs de dispersion (comme par exemple l’écart interquantile).
# Visualiser les premières lignes des données lissées
head(dfBaseLisse)
## Simple feature collection with 6 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 650800 ymin: 6854800 xmax: 652000 ymax: 6855000
## Projected CRS: RGF93 / Lambert-93
## nbObs valeur_fonciere_01 valeur_fonciere_05 valeur_fonciere_09 x y
## 1 55 152000 335000 700000 650900 6854900
## 2 56 152000 319000 675000 651100 6854900
## 3 58 152000 319000 655000 651300 6854900
## 4 60 152000 319000 655000 651500 6854900
## 5 58 182000 324000 655000 651700 6854900
## 6 60 185000 324000 655000 651900 6854900
## geometry
## 1 POLYGON ((650800 6855000, 6...
## 2 POLYGON ((651000 6855000, 6...
## 3 POLYGON ((651200 6855000, 6...
## 4 POLYGON ((651400 6855000, 6...
## 5 POLYGON ((651600 6855000, 6...
## 6 POLYGON ((651800 6855000, 6...
# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="valeur_fonciere_01",
breaks = "quantile",
nbreaks = 5,
lwd=1,
leg_val_rnd = 0,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Premier décile des prix de vente (rayon de 1500m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
type = "choro",
var="valeur_fonciere_05",
breaks = "quantile",
nbreaks = 5,
lwd=1,
leg_val_rnd = 0,
add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
lwd=4,
col="black",add = TRUE)
mf_layout(title = "Médiane des prix de vente (rayon de 1500m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")
Pour que le quantile lissé ait du sens, il faut qu’un nombre conséquent d’observations ait participé à sa construction. Ainsi, il est conseillé d’utiliser un rayon de lissage suffisamment élevé.
summary(dfBaseLisse_paris$nbObs)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 21 715 1029 1096 1430 2907
# Part de centroïdes dans le quantile lissé a été calculé avec moins de 10 observations
dfBaseLisse_paris %>% st_drop_geometry() %>% group_by(nbObs<10) %>% count()
## # A tibble: 1 × 2
## # Groups: nbObs < 10 [1]
## `nbObs < 10` n
## <lgl> <int>
## 1 FALSE 2672
Dans cette section nous allons produire des indicateurs sur des zones à façon.
Pour commencer, nous calculerons le nombre et le prix moyen des transactions immobilières en 2021 dans le “triangle d’or” de la capitale. Ensuite, nous mobiliserons les données carroyées à 200 mètres disponibles sur insee.fr) pour calculer la proportion de logements sociaux autour du White.
Remarque : une fonction généralisant ce second exercice est disponible en ligne, sur insee.fr, en accompagnement des données carroyées.
La zone à façon que nous considérons est le cercle de 1 km de rayon autour du White, à Montrouge. L’objectif est de déterminer la part de logements sociaux dans cette zone, à partir des données carroyées (carreaux de 200 mètres).
Dans un premier temps, on importe la zone à façon. Ici, il s’agit du quartier appelé «~Triangle d’or~», situé dans le 8e arrondissement de Paris. Ce quartier se caractérise par ses immeubles cossus, ses hôtels et ses commerces de luxe.
object <- "diffusion/projet_formation/r_lissage_spatial/triangle_or.kml"
triangleOr <- st_read(paste0(url_bucket,bucket,"/",object))
## Reading layer `triangle_or' from data source
## `https://minio.lab.sspcloud.fr/kantunez/diffusion/projet_formation/r_lissage_spatial/triangle_or.kml'
## using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 2.300717 ymin: 48.86447 xmax: 2.310474 ymax: 48.87203
## Geodetic CRS: WGS 84
leaflet(data=triangleOr) %>% addTiles() %>% addPolygons()
Remarque : ce fichier géogrpahique a pour extension .kml. Ce format se lit sans problème avec la fonction sf::st_read. Pour information, ce triangle a été crée manuellement sur le site du Géoportail.
Pour connaître le nombre et le prix moyen des transactions immobilières réalisées en2021 dans ce quartier, nous procédons par intersection géographique entre des points (les transactions) et un polygone (le quartier). Attention à bien s’assurer que les deux couches géographiques soient bien projetées dans le même système de projection.
# Système de projection de triangleOr
st_crs(triangleOr)[["input"]]
## [1] "WGS 84"
# Transformation du système de projection de triangleOr
triangleOr <- st_transform(triangleOr,crs = 2154)
# Intersection géographique avec les transactions
pos_intersec <- st_intersects(triangleOr,sfBase) %>% unlist()
transac_triangle <- sfBase[pos_intersec,]
# Visualisation de points retenus
transac_triangle
## Simple feature collection with 20 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 648767.5 ymin: 6863055 xmax: 649194 ymax: 6863601
## Projected CRS: RGF93 / Lambert-93
## First 10 features:
## id_mutation date_mutation type_local nombre_pieces_principales valeur_fonciere
## 18629 2021-489256 2021-01-15 Appartement 5 3300000
## 18676 2021-489336 2021-01-25 Appartement 3 590000
## 18677 2021-489337 2021-01-28 Appartement 4 1511345
## 18698 2021-489369 2021-02-05 Appartement 6 2450000
## 18723 2021-489406 2021-01-20 Appartement 3 997000
## 18725 2021-489410 2021-02-12 Appartement 3 1605000
## 18748 2021-489452 2021-02-24 Appartement 2 821500
## 18757 2021-489464 2021-02-04 Appartement 2 1200000
## 18768 2021-489483 2021-02-19 Appartement 2 74900
## 18777 2021-489495 2021-02-19 Appartement 4 2500000
## surface_reelle_bati geometry
## 18629 180 POINT (648918.9 6863145)
## 18676 53 POINT (648903.5 6863601)
## 18677 65 POINT (649194 6863377)
## 18698 133 POINT (648814.5 6863055)
## 18723 72 POINT (648870.5 6863135)
## 18725 100 POINT (648767.5 6863434)
## 18748 60 POINT (648771.3 6863503)
## 18757 74 POINT (648767.5 6863434)
## 18768 44 POINT (648814.2 6863144)
## 18777 100 POINT (648970.8 6863181)
# Cartographie des points dans le triangle
leaflet() %>%
addTiles() %>%
addPolygons(data = triangleOr %>% st_transform(4326)) %>%
addMarkers(data = transac_triangle %>% st_transform(4326))
Remarque 1 : pour tracer cette dernière carte, nous devons re-projeter le triangle et les points en WGS84 (4326), car la fonction Leaflet ne le fait pas automatiquement (contrairement à mapview).
Remarque 2 : on ne voit que 12 points sur la cartes, alors qu’on a trouvé 20 transactions dans le triangle : ceci s’explique car certaines transactions sont géolocalisées au même endroit (appartements situés dans un même immeuble par exemple).
L’ensemble des transactions ayant été récupérées par intersections géographiques, nous pouvons enfin calculer les indicateurs souhaités.
# Nombre de transactions en 2021 :
nrow(transac_triangle)
## [1] 20
# Prix moyen
mean(transac_triangle$valeur_fonciere)
## [1] 1637256
Le principe de la démarche consiste à : 1. Déterminer, pour une zone donnée, l’ensemble des carreaux qui la recouvrent 2. Calculer les agrégats sur cet ensemble de carreaux.
Remarque : L’ensemble de carreaux étant en général plus large que la zone, les agrégats obtenus seront des estimations qui tendent à surestimer les valeurs réelles. Le problème ne se pose pas si vous travaillez directement sur des données disponibles au niveau “individu statistique” (“avec des x,y”).
@Julien : on peut aussi donner un exemple d’agrégation aux x,y avec les données de DVF juste avant l’exemple au carreau. Ca me paraitrait intéressant.
# library(stringr, warn.conflicts = FALSE)
# Charger les données carroyées de Filosofi 2015 à 200m
# Uniquement en région parisienne
# En ne conservant que 6 variables dont le nombre de logements sociaux
st_read_maison <- function(chemin_tab){
requete <- "SELECT IdINSPIRE,Depcom,I_est_cr,Men, Log_soc, geom
FROM Filosofi2015_carreaux_200m_metropole
WHERE SUBSTR(Depcom, 1, 2) IN ('75','92','93','94') "
sf::st_read(chemin_tab, query = requete)
}
object = "diffusion/projet_formation/r_lissage_spatial/Filosofi2015_carreaux_200m_metropole.gpkg"
carreaux <- st_read_maison(paste0(url_bucket,bucket,"/",object))
head(carreaux)
## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 652623.2 ymin: 6858353 xmax: 655268.9 ymax: 6866380
## Projected CRS: RGF93 / Lambert-93
## IdINSPIRE Depcom I_est_cr Men Log_soc
## 1 CRS3035RES200mN2893400E3763200 75119 0 990 937
## 2 CRS3035RES200mN2890000E3762200 75111 0 926 80
## 3 CRS3035RES200mN2893400E3762400 75119 0 508 323
## 4 CRS3035RES200mN2891800E3763600 75119 0 633 0
## 5 CRS3035RES200mN2890800E3763200 75120 0 349 188
## 6 CRS3035RES200mN2885800E3760600 75113 0 751 344
## geom
## 1 MULTIPOLYGON (((654521.9 68...
## 2 MULTIPOLYGON (((653843.4 68...
## 3 MULTIPOLYGON (((653725.1 68...
## 4 MULTIPOLYGON (((655069.7 68...
## 5 MULTIPOLYGON (((654764.7 68...
## 6 MULTIPOLYGON (((652641.8 68...
# Création du cercle de 1km autour du White
white <- st_sfc(st_point(c(649218.36,6857569.14)),crs=2154)
buffer_white <- st_buffer(white,1000)
# Sélection des carreaux intersectant notre zone d'intérêt
carreaux_select <- carreaux[unlist(st_intersects(buffer_white,carreaux)),]
# Cartographie
leaflet() %>%
addTiles() %>%
addPolygons(data=carreaux_select %>% st_transform(4326),
label = "Carreaux sélectionnés") %>%
addPolygons(data=buffer_white %>% st_transform(4326) ,
color="red",
fill=F,
label = "Zone à façon")
Cette carte représente les carreaux issus de Filosofi carroyé en 2015 qui intersectent notre zone d’intérêt. À noter que ces carreaux sont « penchés » : ceci est la résultant de leur reprojection en Lambert 93. En effet, ces carreaux ont été produits en projection LAEA (standard européen), ce qui leur donne cet aspect une fois reprojetés dans le standard français.
Maintenant que nous avons sélectionné ces carreaux, il suffit de calculer l’indicateur souhaité.
tx_logSoc_white <- sum(carreaux_select$Log_soc)/sum(carreaux_select$Men)
cat("On compte ",
round(tx_logSoc_white*100),
"% de logements sociaux dans un rayon de 1km autour du white")
## On compte 33 % de logements sociaux dans un rayon de 1km autour du white
Vous êtes désormais capables de carroyer et lisser toutes les données que vous aurez à votre disposition !
Reproducibilité
sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] leaflet_2.1.0 aws.s3_0.3.21 mapview_2.10.0 mapsf_0.4.0
## [5] btb_0.1.30.3 sf_1.0-6 dplyr_1.0.8 unilur_0.4.0.9100
## [9] knitr_1.37
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.2 jsonlite_1.8.0 RcppParallel_5.1.5
## [4] assertthat_0.2.1 sp_1.4-6 highr_0.9
## [7] stats4_4.1.2 renv_0.15.2 yaml_2.3.5
## [10] pillar_1.7.0 lattice_0.20-45 glue_1.6.2
## [13] uuid_1.0-3 digest_0.6.29 RColorBrewer_1.1-2
## [16] colorspace_2.0-3 leaflet.providers_1.9.0 htmltools_0.5.2
## [19] pkgconfig_2.0.3 raster_3.5-15 purrr_0.3.4
## [22] scales_1.1.1 webshot_0.5.2 brew_1.0-7
## [25] svglite_2.1.0 terra_1.5-21 satellite_1.0.4
## [28] tibble_3.1.6 proxy_0.4-26 generics_0.1.2
## [31] farver_2.1.0 ellipsis_0.3.2 cli_3.2.0
## [34] magrittr_2.0.2 crayon_1.5.0 evaluate_0.15
## [37] fansi_1.0.2 xml2_1.3.3 class_7.3-19
## [40] tools_4.1.2 lifecycle_1.0.1 stringr_1.4.0
## [43] munsell_0.5.0 compiler_4.1.2 jquerylib_0.1.4
## [46] e1071_1.7-9 systemfonts_1.0.4 rlang_1.0.1
## [49] classInt_0.4-3 units_0.8-0 grid_4.1.2
## [52] leafpop_0.1.0 htmlwidgets_1.5.4 aws.signature_0.6.0
## [55] crosstalk_1.2.0 leafem_0.1.6 base64enc_0.1-3
## [58] rmarkdown_2.11 codetools_0.2-18 DBI_1.1.2
## [61] curl_4.3.2 R6_2.5.1 fastmap_1.1.0
## [64] utf8_1.2.2 KernSmooth_2.23-20 stringi_1.7.6
## [67] Rcpp_1.0.8 vctrs_0.3.8 png_0.1-7
## [70] tidyselect_1.1.2 xfun_0.29